function [A, Z, objHistory] = JMKSC(H,k,alpha,gamma,beta)
[~,n,nKernel]=size(H);
sumH=zeros(n);
for i=1:nKernel
    sumH=sumH+H(:,:,i);
end
commonH=sumH/nKernel;
tol = 1e-2;
maxIter = 100;
one = ones(n,1);    %nx1
I = eye(n);         %nxn
Z=eye(n);
S = zeros(n);
A = Z;
L = diag(A*one)-A;
w = zeros(nKernel,1);
for j=1:nKernel
    w(j)= correntropy(H(:,:,j),commonH);
end
iter = 0;
alphaovergamma = alpha/gamma;
objHistory = [];
objHistory = [objHistory; CacaulateObj(H, commonH, Z, L, S, A, w, alpha, beta, gamma)];
while iter < maxIter
    iter = iter + 1;
    %********************************************
    % update Z
    %********************************************
    Zk = Z;
    Z = (gamma*I + commonH)\(commonH+gamma*A);
    
    %********************************************
    % update S
    %********************************************
    [V, D] = eig(L);
    D = diag(D);
    [~, ind] = sort(D);
    S = V(:,ind(1:k))*V(:,ind(1:k))';

    
    %********************************************
    % update A
    %********************************************
    Ak = A;
    A = Z-alphaovergamma*(repmat(diag(S),1,n)-S);
    A = max(0,(A+A')/2);
    A = A-diag(diag(A));
    L = diag(A*one)-A;
    
    %********************************************
    %update commonH and w
    %********************************************
    commonHk=commonH;
    sumH=zeros(n);
    sum1=0;
    w = zeros(nKernel,1);
    for j=1:nKernel
        w(j)= correntropy(H(:,:,j),commonH);
        sumH=sumH+w(j)*H(:,:,j);
        sum1=sum1+w(j);
    end
    commonH=(2*Z'-Z*Z'-eye(n)+2*beta*sumH)/(2*beta*sum1);
    commonH = commonH .* (commonH > 0);
    commonH=(commonH'+commonH)/2;
    
    
    diffZ = norm(Z-Zk,'fro');
    diffA = norm(A-Ak,'fro');
    diffH = norm(commonHk-commonH,'fro');
    objHistory = [objHistory; CacaulateObj(H, commonH, Z, L, S, A, w, alpha, beta, gamma)];
    stopC = max([diffZ,diffA,diffH, objHistory(end-1) - objHistory(end)]);
    if (iter>5) && (stopC <  tol)
        break
    end
end
end


function obj = CacaulateObj(H, commonH, Z, L, S, A, w, alpha, beta, gamma)
f1 = trace(commonH - 2*commonH*Z+Z'*commonH*Z);
f2 = sum(sum(S.*L));
f3 = 0;
for j=1:size(H,3)
    f3 = f3 +w(j)*norm(H(:,:,j)-commonH,'fro')^2;
end
f4 = .5 * norm(Z-A,'fro')^2;
obj = 0.5*f1 + alpha * f2 + beta * f3 + gamma * f4;
end

function value = correntropy(H1,H2)
Ha = H1-H2;
[m1,n1]=size(Ha);
d2 =  (norm(Ha,'fro')^2)/(m1*n1);
Ha = Ha.*Ha;
Ha = exp(-(Ha)/(d2));
value= sum(Ha(:))/(m1*n1);
value = 1-sqrt(1-value);
end